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Abstract 

We show that the standard class of optimal control models in OCMat 1 can be used to analyze 
ID spatial distributed systems. This approach is an intermediate step on the way to the FEM 
discretization approach presented in trass and Uccker [2015]. Therefore, the spatial distributed 
model is transformed into a standard model by a finite difference discretization. This (high dimen¬ 
sional) standard model is then analyzed using OCMat. As an example we apply this method to the 
spatial distributed shallow lake model formulated in Brock and Xepapadeas [2008]. The results are 
then compared with those of the FEM discretization in Grass and Uecker [2015]. 

Keywords: spatial distributed optimal control model, finite difference discretization, shallow lake model, 
patterned indifference threshold point 


1 Introduction 

The analysis of spatial distributed optimal control problems (over an infinite time horizon) become 
an important issue in economic modeling, see e.g., Brock and Xepapadeas [2008], Brock et al. [2014]. 
In technical terms this means that the time evolution of the space distributed states are described by 
parabolic PDEs. In 3rock and Xepapadeas [2008] the authors provide a local stability analysis of the 
equilibria, i.e. the solutions of the elliptic PDEs associated to the derived canonical system. They 
derive a set of conditions that cause a Turing instability of these equilibria and call the bifurcating 
patterned equilibria POSS solutions, i.e. Patterned/Heterogeneous Optimal Steady States in contrast 
to Homogeneous or Flat Optimal Steady States (FOSS). 

But, as is shown in Grass and Uecker [2015], a local stability analysis is not sufficient to prove the 
optimality of heterogeneous equilibria solutions. To answer the question of optimality the objective 
values of the paths that converge to the different equilibria have to be compared. This analysis 
also sheds new light on the discussion about indifference threshold and threshold points, cf. Kiseleva 
and Wagener [2010], Kiseleva [201 ]. We therefore introduce a new terminology of defective and 
non-defective equilibria. This properties distinguishes between optimal equilibria that are stable or 
unstable. 

Since in general the PDEs cannot be solved analytically we have to resort to numerical methods. 
In frass and Uecker [2015] we present a numerical procedure relying on a FEM discretization of the 
derived canonical system combined with a continuation strategy, analogous to the approach in OCMat, 
cf. Grass [2012]. As an example we used the distributed shallow lake model Eq. (19). 

1 OCMat is a MATLAB package that provides tools for the numerical analysis of (non-distributed) optimal control 
problems, specifically over an infinite time horizon. It can be downloaded from http://orcos.tuwien.ac.at/research/ 
ocmat_software. 
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This allowed us to identify parameter regions where non-defective POSS exist, and where defec¬ 
tive POSS and the according stable manifolds separates the regions of attractions of FOSS and POSS 
(threshold points/manifolds). Additionally we were able to show the existence of homogeneous/hetero¬ 
geneous indifference threshold (Skiba) points and calculated the connecting manifold of indifference 
threshold points. 

In this paper we take an intermediate step to the full FEM approximation of the canonical system, 
to demonstrate the capabilities of the standard optimal control class of OCMat, i.e. infinite time horizon 
problems, where the time evolution of the states is described by ODEs. To apply the OCMat processes 
for the initialization and file generation directly we discretize the PDEs of the state equations by a 
finite difference scheme (FD). This yields a number of ODEs that can be handled by OCMat. The 
numerical results of this approach are then compared with the analogous results of trass and Uecker 
[2015]. 

The article is structured in the following way. We start with a discussion and introduction of general 
terms and shortly present the finite difference discretization in Section 2. In Section 3 we summarize 
important properties and results of the 0D (non-distributed) shallow lake model, introduced and 
analyzed in Scheffer [1998], Maler et al. [2003], Carpenter and Brock. [2004] and Wagener [2003]. 
In the next Section 4 the ID spatial distributed shallow lake model is formulated together with its 
discretized counterpart. The latter model is then numerically analyzed in detail. 

2 General Definitions 


2.1 Models of spatial dimension 0 (OD model) 


r oo 

max / e~ pt g(x(t), u(t))dt 
«(■) Jo 

(la) 

s.t. x{t) = f(x(t),u{t)) 

(lb) 

x(0) = xq e M n . 

(lc) 

with / € CPfW 1 x g € C 2 (M n x 

Let (x*(-), u*(-)) be an optimal solution of Eq. (1). Then there exists A(-) such that 
is a solution of the canonical system 

,«*(•), A(-)) 

x(t) = f(x*{t),u{t)) 

(2a) 

A(t) = p\(t) -'H(x*(t),u*(t),\(t),\ 0 ) 

(2b) 

x*{0) = x 0 

(2c) 

with 


u*(t) = argmax%(x*(f), u, A(f), Ao) 

U 

(2d) 

and 


H{x, u, A, A 0 ) := A 0 g(x, u ) + A T /(x, u), 

(2e) 


To ease the notation we make the following assumptions 

Assumption 2.1. Problem Eq. (1) is normal, i.e., Ao = 1. Therefore, we omit the argument Ao- 

Assumption 2.2. Let (x*(-), «*(•)) be an optimal solution and A(-) the according costate. Then, there 
exists an explicit function 

u°(x, A) € C 2 (IT x 

such that for every t 

H(x*(t ), u°(x*(t), A (t)), A(t)) = ma x'H(x*(t),u, A (t)). 

U 
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Then an optimal solution (x*(•), u*(•)) can be found among the solutions satisfying 

x(t) = f(x(t),X(t)) 

A (t) = pX(t) - H x (x(t ), A (t)) 
x(0) = x 0 . 

with 

:= f(x(t),u°(x(t),X(t))) 

and 

TL{x{t), X(t)) := TL{x{t), u°(x(t), X(t)), X(t), 1). 


(3a) 

(3b) 

(3c) 


(3d) 


Subsequently we will omit the bar sign. 


Definition 2.1 (OSS). Let (x*(■), u*(■)) with x*(-) = x E and u*(-) = u 6 M m be an optimal 
solution of problem (1) with x(0) = x . Then (x,u) is called an optimal equilibrium and denoted as 
OSS. 


Let (x, A) E 
as CSS and 


be an equilibrium of the canonical system Eq. (3). Then (x, A) E M 2n is denoted 


J(x, A) : = 


t 

d/(x, A) 

df(x, A) \ 


dx 

dA 

V 

dH°{x,\) 

dU°(x,\) | 

dx 

dA / 


( 4 ) 


(x,\) 


is the according Jacobian matrix, and if there is no ambiguity we simply write J. The eigenspaces 
corresponding to J(x, A) are denoted as 


E s (i,A) :={(eC: J{x, X)v = £v with Re£ < 0}, n s := dimE s (x, A) (5a) 

E u (f, A) := {£ E C : J(x, X)v = £,v with Re£ > 0}, n u := dimE u (x, A) (5b) 

E c (x, A) := {£ E C : J(x,A)u = 0}, n c := dimE c (x, A). (5c) 

Definition 2.2 (Saddle Point Property). Let ( x, X ) E M 2n be an equilibrium of Eq. (3). If 

dimE s (x, X) = n (6) 


then it is said, that the equilibrium satisfies the saddle point property (SPP). The equilibrium (x, A) 
is denoted as CSSP. Otherwise it is denoted as CSS~. The number 


d(x, A ):=n s -n u -n c (7) 

is called the defect of (x, A). An equilibrium with defect d(x, A) < 0 is called defective, otherwise it is 
called non-defective. If (x, A) is defective and the according (x,u°(x, X)) is OSS, then (x, u°(x, A)) is 
called defective otherwise it is called non-defective. 

Proposition 2.1. Let (x, X) E M 2n be an equilibrium of Eq. (3) and p> 0. Then (x, X) satisfies the 
saddle point property iff every eigenvalue £ of the according Jacobian J(x, A) satisfies 

|Ree - f I > f (8) 

Proof. In Grass et al. [2008] it is proved that there exist n (not necessarily distinct) complex numbers 
£ E C with Re£ > 0 such that any eigenvalue £ of the according Jacobian, satisfies 

H+ f " or «=H 
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This symmetry together with SPP yields 

|IteC-fl = Ree> f 

and the last inequality is identical to 

Re£>| iff Re£ < 0. 

Therefore, n eigenvalues have a negative real part finishing the proof. □ 

Remark 2.1. Proposition 2.1 allows us to formulate the SPP by the equivalent Eq. (8). A definition of 
SPP claiming Eq. (8) has the advantage that it also can be applied to equilibria of distributed systems. 
Where a definition relying on the dimension of the stable manifold fails. 


2.2 Models of spatial dimension 1 (ID model) 


We assume that the spatially distributed model is derived by the introduction of a diffusion term, 
whereas the functions / and g are the same as for model (1). This yields 


r*oo rL 


max 




“(v) Jo J-L 


g(x(z, t),u(z , t))dzdt 


s.t. = f(x(z,t),u(z,t)) + D d 


dx(z , t ) 


dz 


= 0 


±L 


x(z, 0) = xo(z), ze[—L,L\. 

Or transforming [-L, L] into [0,1] yields 

roo rl 

max / / e~ pt g(x(z,t),u(z,t))dzdt 

“(v) Jo Jo 

d . ....... D d 2 x(z,t ) 

s.t. — x(z, t) = f{x{z, t),u(z, t)) + dz2 


dx(z,t) 


dz 


= 0 


0,1 


x(z,0) = x 0 (z), z€[0, 1]. 


(9a) 

(9b) 

(9c) 

(9d) 


Applying Pontryagin’s Maximum Principle for PDEs, see e.g., Troltzsch [2009], we can derive, analo¬ 
gous to Eq. (3), the canonical system for (9) as 


d_ 

dt 

d_ 

dt 


t) = f(x(z, t), A (z, t)) + D 


d 2 x(z, t) 
dx 2 


^ dU{x(zfi),\(zfi)) 
\{z,t) = pX{z,t )- 7 ^- D 


d 2 X(z,t ) 
dx 2 


d n x(z,t) | 0)1 = 0 
d n X{z,t)\ 01 = 0 


x(z,0) = x 0 (z), z G [0,1]. 


(10a) 

(10b) 

(10c) 

(lOd) 

(10e) 


For the numerical analysis we can than e.g. use a finite element method (FEM) for the discretization 
of Eq. (10). For the distributed shallow lake model (cf. Section 4) this has been carried out in Grass 
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and Uecker [2015]. In this article we want to demonstrate the capabilities of OCMat and transform it 
into a high dimensional 0D model. 

For the FDM discretization we use an equidistant grid of length h with Nh = 1 and N E N. 


d (^ 
- X (z) 


Zi 


dz 2 

f 

Jo 


x(z) 


x(zi + h) — x(zi — h) 

2 h 

x{zi — h) - 2x(zi) + x(zi + h ) 


Zi 


h 2 


SN-l 




(11a) 

(lib) 

(11c) 


Thus, model (9) is approximated by 


s.t. 


with 


1 

max — 
u 0 {-),-,u N (-) N 


fN-l 


/ 0 


Pt 9{xi{t),u,i{t)) + 


g(x 0 (t),u 0 (t)) + g(x N (t),u N (t)) 


I d t 


\ 1=1 
DN 2 

Xi(t) = f(xi(t),Ui(t )) + 9 - 2 Xi(t) +x i+ i(t)) 

(2 L) 2 

X\ (t) - x-i(t) = x N+ i(t) - XN-I(t) = 0, t > 0 
Xi(0) = x ifi . 


I 

Zi= N' 7 = 0 ’---’ iV 

Xi(t) := x(zi,t), Ui(t) := u(zi,t) 

Xi ■■= x(Zi,-), Ui'.= u(Zi,-) 


(12a) 

(12b) 

(12c) 

(12d) 


Note that we do not consider the problem under which conditions model (12) approximates (9). We 
rather take, in the specific case of the spatial shallow lake model, model (12) for granted to see if 
OCMat can handle such a problem and compare the results with those from Grass and Uecker [2015]. 
The canonical system for model (12) becomes 


and 


Xi(t) = f(xi(t),Xi(t)) + V[ x \t) 

A i{t) = p\i(t) - n x (xi(t), A i(t)) - X>- A) (f) 

Xi( 0) = Xi i0 . 


D 


Ui(xi , A i) 

DN 2 

(2 W 




2D(x\ - x Q ) 

D{xi -1 - 2 Xi + x i+ i) 
2D(x N -i - x N ) 

f D( Ai-2A 0 ) 

D{ 2Ao — 2Ai + A 2 ) 

D(Xi-i - 

D(Xn-2 
D(Xn-i 


i = 0 

i = 1,..., N — 1 
i = N 

i = 0 
i = 1 

i = 2,.. . ,N — 2 


2A % + Aj+i) 

- 2Ajv —1 + 2A n) i = N — 1 

- 2A n i = XI) 


(13a) 

(13b) 

(13c) 
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To abbreviate notation we introduce 


x d := (xq , • • •, xJf) T G ]^ n ( Ar + 1 ) 

x d -.= (\t,...,\ t n ) t €r ( Ar+1 ) 

u d -.= ( ? 4,...,4) T Grw. 

Definition 2.3 (FOSS and POSS). Let (x d, *(-), it d, *(-)) with x d, *(-) = x d € M n ( 7V+1 ) and u d, *(-) = 
u d G M m ( Ar+1 ) be an optimal solution of problem (1) with x d (0) = x d . If 

Xq = X\ = ■ ■ ■ = X TV = X G M n 

then ( x d ,u d ) is called a flat (homogeneous) optimal steady state (FOSS), otherwise it is called an 
patterned (heterogeneous) optimal steady state (POSS). 

Definition 2.4 (FCSS and PCSS). Let ( x d ,X d ) G M 2n ( ,v + 1 ) be an equilibrium of the canonical system 
Eq. (13). Then (x d ,X d ) is called a flat (homogeneous) steady state (FCSS) iff 

xq = x\ = ■ ■ ■ = xn = x (14) 


otherwise it is called a patterned (heterogeneous) steady state (PCSS). If a FCSS (PCSS) satisfies 
SPP it is denoted as FCSSP (PCSSP) otherwise it is denoted as FCSS~ (PCSS~). 

Definition 2.5 (State-Costate space). Let (x d (■), \ d (■)) be a solution of the canonical system Eq. (13). 
Then the representation (||x d ||(-), ||A d ||(-)) with 




llv?(*)ll + 


mm + Km 


y = x, or y = A, j = 1,... ,n. 


(15) 


is called a solution path in the state-costate space. 

Remark 2.2. In Brock and Xepapadeas [2008] FOSS (POSS) were introduced as flat (patterned) 
optimal equilibria of the canonical system satisfying SPP. We enhanced this terminology for two 
reasons 


1. For better clearness we decided to make a further distinction between the canonical and optimal 
system. 

2. There exist optimal equilibria that do not satisfy SPP, cf. Section f.f.l. 


3 The shallow lake model without spatial diffusion 

A well known version of the shallow lake model, see e.g. Wagener [2003], can be formulated as 

roo 

max / eT pt (ln(w(£)) — cP(t) 2 )dt (16a) 

«(•) Jo 

s.t. P(t) = u(t) - bP(t) + i (16b) 

P(0) = P 0 > 0. (16c) 

We sometimes refer to model (16) as 0D shallow lake model. 
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By Pontryagin’s Maximum Principle we find the canonical system 


with 


m = u°(t) - bp(t) + jTLL 

(17a) 

A(f) - 2cP{t) + \{t) (p + b (1 + p ^ 2 ) 

(17b) 

(0 

II 

(17c) 

(V 

(17d) 


Let (P*(■), u*(•)) be the optimal solution of Eq. (16); then P*(-) is the unique solution of the IVP 

Pit? 


P(t) = u*(t) — bP[t ) + 
P( 0) = P Q . 


1 + P{t? 


(18a) 

(18b) 


The ODE (18a) is called the optimal system for u*(t). In Vagener [2003] it is proved that every 
optimal solution (P*(■), u*(•)) for arbitrary Pq > 0 converges to an equilibrium ( P,u ) of the optimal 
system with u > 0, usually depending on Po. 

A detailed bifurcation analysis in the parameter space (6,c), see, e.g., Wagener [2003], reveals the 
existence of regions in the parameter space where the optimal system consists of 

• One globally stable optimal equilibrium ( P,u ). 


• Two locally stable equilibria ( P 0 ,u 0 ) and ( P e ,u e )• These are separated by one of the following 
state values 


— The state value P u of an unstable optimal equilibrium ( P u ,u u )• 

— An indifference threshold point Pj also called Skiba point. 

Thus, we can give a full classification of the optimal solutions. In the case that three optimal equilibria 
exist we choose P Q and P e such that P Q < P u < P e . There are intermediate cases (bifurcation cases) 
where equality holds that are not specifically mentioned. We also refer to P e as the eutrophic and to 
P a as the oligotrophic equilibrium. 

Global stable: For any initial state Po > 0 there exists a unique solution (P*(-), u*(-)) that converges 
to ( P,u ), which is independent of Po- See Fig. 2c. 

Local stable I: For any initial state 0 < Po < P u there exists a unique solution (P*(-), «*(•)) that 
converges to (P 0 ,ui), which is independent of Po. For Po > P u there exists a unique solution 
(P*(-),«*(•)) that converges to (P e ,U2), which is independent of Po- For Po = P u the optimal 
solution is (P*(■), u*(■)) = ( P u ,u u ). See Fig. 3c. P u is called a threshold point. 

Local stable II: The first two statements of the previous case remain true, replacing P u by Pj. For 
Po = Pi there exist two optimal solutions (P*(•), u* (•)), i = 1,2 that converge to (P 0 ,ui) and 
(P e , uf) with P Q < Pj < P e . See Fig. 4c. Pi is called an indifference threshold point or Skiba 
point. 

From the perspective of optimal control theory the last two cases are of specific interest. The first 
case is often referred to as history dependence, i.e., the optimal solution depends on its initial starting 
point. In the second case we additionally observe the non-uniqueness of the optimal solution, i.e., 
indifference. For a detailed discussion and description of the underlying bifurcations we refer to 

Kiseleva and Wagener [2010], Kiseleva [2011]. 

The parameter values for these two prototypical cases are specified in Table 1. In the first case we 
find an indifference threshold point and in the second case a threshold point. 
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Model (16) 


Model (19) specific 


Scenario 

P 

c 

b 

D 

L 

N 

I 

0.03 

0.5 

0.65* 

0.5 

2tt/0.44 

50 

II 

0.3 

3.5* 

0.55 

0.5 

2tt/0.44 

50 


Table 1: The parameter values for the two considered scenarios for the OD and ID model. The values 
with the superscript * denote the free parameter. 


Bifurcation-analysis Anyhow, this classification is the result of an intensive numerical analysis of 
the canonical system Eq. (17). This analysis covers a bifurcation analysis of its equilibria (see Fig. 1) 
and the calculation of the related stable paths and their objective value. The numerical computation 
is necessary since the local properties of an equilibrium (P, A) of Eq. (17) let us not deduce that 
the corresponding equilibrium ( P,u ) with u = 1/A is an optimal equilibrium of the optimal system 
Eq. (18a). This is specifically true in the case that multiple equilibria of the canonical system exist. 
Therefore, there is no one-to-one correspondence between the bifurcations of the optimal and canonical 
system. 

Unique optimal equilibrium The importance of a numerical analysis of the stable paths comes 
specifically clear in scenario I with b = 0.75. In Fig. 2 we see that there exist three equilibria in the 
canonical system (two saddles and one unstable focus), whereas the optimal system only consists of 
one globally stable equilibrium. 

Indifference threshold point For b = 0.65 the number and properties of the equilibria of the 
canonical system remain the same, but the optimal system consists of two locally stable equilibria 
separated by an indifference threshold point P/, cf. Fig. 3. Thus, a local stability analysis of the 
equilibria has to be supported by the global analysis the according stable manifolds. Even if the first 
task can be realized analytically, the stable paths can only be calculated analytically in very rare cases. 
Usually and specifically in our case we have to resort to numerical methods to solve the latter task. 

Indifference point For scenario II with c = 3.5 the canonical system exhibits two saddles and one 
unstable node (cf. Fig. 4a). Calculating the stable paths and comparing the objective values we find 
that the unstable equilibrium is a threshold point (cf. Fig. 4b). This unstable equilibrium separates the 
regions of attraction for the two locally stable equilibria corresponding to the two saddles(cf. Fig. 4c). 
The second case with c = 3.0825 yields qualitatively the same result and is not depicted. 


4 The shallow lake model with spatial diffusion 


An extension of the shallow lake model (16) to the class of spatially distributed models, proposed in 

Brock and Xepapadeas [2008], is given by 


max / e pt \n(u(x,t)) — cP(x,t)dx dt (19a) 

«(•>•) Jo Jn 

s.t. ^P(x,t) = ii(i,f)-6P(i,f) + ^^ + D-^P(x,t) (19b) 

d n P(x,t) | 9n = 0 (19c) 

P(x,t)\ t=0 = Pq(x), (19d) 


Contrary to 3rock and Xepapadeas [2008] we formulated the problem for Neumann conditions Eq. (19c), 
the so called zero flux boundary condition, instead of the periodic boundary conditions. From an in- 
terpretational point of view we assume the lakes to be located consecutively in a row and there is 
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(a) Equilibria bifurcation diagram (b) Equilibria bifurcation diagram 

Figure 1: The • denote equilibria satisfying the SPP and o denote equilibria not satisfying SPP. 
The figures in the first row show the bifurcation parameter versus the (absolute) state value. In 
the second row the norm of the equilibrium is plotted versus the bifurcation parameter. Panel (a) 
(p = 0.03, c = 0.5 and varying b ) show the existence of two separated branches of equilibria. In the 
interval [0,0.727] there exist three equilibria. Panel (b) (p = 0.3,6 = 0.55 and varying c) shows the 
existence of one branch of equilibria. In the interval [2.566,3.556] there exist three equilibria. 




(a) State-Costate (b) State-Objectivevalue (c) Optimal System 


Figure 2: For p = 0.03, c = 0.5 and 6 = 0.75 there exist three equilibria in the canonical system (a). 
The optimal system (c) only consist of one globally optimal equilibrium. In panel (a) the • denote 
saddles of the canonical system, the o an unstable focus. In panel (c) the • denotes the globally stable 
equilibrium. 
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(a) State-Costate (b) State-Objectivevalue (c) Optimal System 


Figure 3: For p = 0.03, c = 0.5 and b = 0.65 there exist three equilibria in the canonical system 
(a). The optimal system (c) consists of two locally optimal equilibria. The basins of attraction are 
separated by an indifference threshold point (Skiba point) Pj. with a discontinuous dynamics at Pj. 
In panel (a) the • denote saddles of the canonical system, the o an unstable focus. In panel (c) the • 
denote the locally stable equilibria. 



(a) State-Costate (b) State-Objectivevalue (c) Optimal System 

Figure 4: For p = 0.3, c = 3.5 and b = 0.55 there exist three equilibria in the canonical system 
(a). The optimal system (c) consists of two locally optimal equilibria. The basins of attraction are 
separated by the optimal unstable equilibrium P u . In panel (a) • denote saddles of the canonical 
system, o an unstable node. In panel (c) • denote the locally stable optimal equilibria and o denotes 
the optimal unstable node. 
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no flow out of lake. Periodic boundary conditions refer to a ring of lakes. Brock and Xcpapadcas 
argue that they use periodic boundary conditions to exclude effects induced by the conditions at the 
end points. Anyhow, since this point does not touch our argument that it is necessary to analyze the 
global behavior of the optimally controlled system, we changed the model formulation in this respect. 

Using the (FDM) discretization proposed in Section 2.2 we find 


with 


max 
u 0 {-),-,u N (-) [Jo 


“OO 

e~ rt G(P 0 (t),..., P N (t), u 0 (t ),.. 

1 

• ,u N (t )) dtj 

(20a) 

p.(f) 2 

bm + TTW)2+ D ( p i-iW- 

2 Pi(t) + Pi+l(t)) 

(20b) 

= PN+i(t) — Pn-i{ t) = 0, t > 0 


(20c) 

3 


(20d) 


x i = i = 0,...,N, D := D 


Pi(t ) := P(xi,t), Ui(t ) := u(xi,t) 

Pi := P{xi,-), Ui:=u(xi,-) 

n , T3 D N 1 ( r> \ , 9(PoiUo) + g{PN,UN) 

G{P 0 , • • • , Pn, U o, • • ■ , UN) := T7 2^ Ui > ^ -9- 


2—1 


g(Pi,Ui ) := ln(rti) - cPf 


Applying Pontryagin’s Maximum Principle on Eqs. (11) yields the canonical system 

\2 


m = u °(t ) - bm) + , + v\ p \t) 


i +m) 2 


Xi(t) = CiPi(t ) + A i(t) ^p + 6- 

Pi(0) = Pi,o > 0, i = 0,...,N 

1 


2 Pi(t) 


a +m?f 


- 


u°(t) = < 


2 A i(t) 


Ci : = 


: = 


K(t) 

c i = 0, N 


2c i = l,...,N -1 


i = 0, N 
i = 1,..., N — 1 


i = 0 


2D(P 1 (t) - P 0 (t)) 

D(Pi-i(t) - 2 Pi(t) + P i+ i(t)) i = 1,..., N - 1 
[2D(P N _ 1 (t) - P N (t)) i = N 


V[ x \t) := 


D(Xi(t) — 2A 0 (t)) 

i = 0 

D{2Xo(t) — 2Ai(f) + X 2 (t)) 

i = 1 

D(Xi~i(t) — 2A i(t) + Xi + i{t)) 

i = 2 ,..., 

D{^N- 2 {t) — 2Aw-i(t) + 2Ajv(i)) 

i = N — 1 

D{^N-i{t) — 2A N{t)) 

i = N 


N -2 


(21a) 

(21b) 

(21c) 

(21d) 


In Appendix B the details for an implementation of the model (20) in OCMat are explained. 
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4.1 Equilibria of the canonical/optimal system 

There is an intimate relation between the CSS of the canonical system Eq. (17) and FCSS of the 
canonical system Eq. (21) 

Corollary 4.1. Let (P d ,u d ) G M 27V+2 be FOSS then {P d ,\/{2uq)) is an equilibrium of the canonical 
system Eq. (17). 

Corollary 4.2. Let (P, X) be an equilibrium of the canonical system Eq. (17). Then P d := (P ,..., P) 
and X d := (2A, A,..., A, 2A) is a FCSS. 

Corollary 4.3. Let (P, A) be a saddle of the canonical system Eq. (17). Then for D small enough 
(. P d , X d ) G M 27V+2 defined in Corollary f.2 is FCS&. 

PCSS can in general not be calculated analytically therefore we have to resort to numerical meth¬ 
ods. Using Corollary 4.2 we can start a bifurcation analysis of Eq. (21) with an equilibrium of Eq. (17). 
The PCSS emerge from branching points of the bifurcation curve. In Grass and Uecker [2015] the ac¬ 
cording bifurcation analysis is done using pde2path, a MATLAB package for the bifurcation analysis 
of elliptic PDEs, see Uecker et al. [2014] and Dohnal et al. [2014]. Since the actual model (16) is a 0D 
optimal control model with a finite number of states (N + 1) the bifurcation analysis of Eq. (21) is 
done using a modified version of CL_MATC0NT 2 . 

We analyzed the two different scenarios specified in Table 1. 

First Scenario The according bifurcation analysis with respect to b is depicted in Fig. 5a. The 
black curves represent the bifurcation curves of the FCSS. These curves exhibit the same shape as the 
corresponding bifurcation curves for the 0D model in Fig. la. For the lower branch we additionally 
find four branching points o, where the branches of the PCSS emanate (red, green, magenta and 
cyan). Along the bifurcation curves of the PCSS we find additional branching points and calculated 
the according bifurcation curves (brown, dark green and orange). Thus, for b = 0.65 we find in total 
two FCSS 0 (correspondig to the oligotrophic and eutrophic equilibrium in the 0D model), one FCSS - , 
thirteen PCSS - and one PCSS 0 . In fact the brown and dark green branch consists of two distinct 
bifurcation curves with equilibria that are spatially symmetric. 

Second Scenario The according bifurcation analysis with respect to c is depicted in Fig. 5b. The 
black bifurcation curve of the FCSS consists of one branch, exhibiting two fold bifurcations and six 
branching points o. These six branching points are connected by three bifurcation curves of the PCSS 
(red, magenta and green). The red and green curve consist of two spatially symmetric branches. From 
the magenta PCSS bifurcation curve two further PCSS branches emanate (brown and orange). 

In Section 4.4 we consider two specific cases for c = 3.5 and c = 3.0825. In the first case there 
exist two FCSS 0 (correspondig to the oligotrophic and eutrophic equilibrium in the 0D model), one 
FCSS - , ten PCSS - and two PCSS 0 . In the latter case there exist two FCSS 0 , one FCSS - and four 
PCSS - . 

4.2 First scenario optimal solutions exemplified on two cases 

Most of the results for the specific cases b = 0.75 and b = 0.65 have already been discussed in trass 
and Uecker [2015]. Therefore we only give a brief summary and concentrate on some aspects that 
were easier to get using OCMat. 

For b = 0.75 there exists a single FOSS, structurally reproducing the result of the according 0D 
model. 

For b = 0.65 the main results are 

2 This modified version is available from the author. 
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(a) First Scenario b, p = 0.03, c = 0.5 (b) Second Scenario c, p = 0.3, b = 0.55 

Figure 5: A bifurcation analysis for the canonical system Eq. (21) is depicted in (a) and (b). The 
symbols: • denote equilibria satisfying the SPP, o equilibria not satisfying the SPP and o branching 
points of the equilibria manifold. The FOSS curves are black and the HOSS curves are colored. Along 
the red HOSS curve two fold bifurcations occur, the HOSS candidates between these fold bifurcations 
satisfy the SPP. 


• None of the PCSS are optimal. 

• The two FCSS 0 are optimal. 

• The according basins of attractions are separated by an indifference threshold manifold. 

• In dr ass and Uecker [2015] we computed a homogeneous and patterned indifference threshold 
point. The homogeneous indifference threshold point coincides with the value of the indifference 
threshold point in model (16). 

Next we explain in detail the calculation of a solution converging to an equilibrium satisfying SPP. 
Subsequently the necessary steps for the detection of an indifference threshold point are explained. 
Details for the calculation in OCMat can be found in Appendix B. 

4.3 A locally optimal patterned equilibrium 

To prove that PCSS 0 (-Pp CSS o, ^pcss°)’ see Fig. 5a, is an optimal equilibrium we have to show that 
there exist no other solution (P d (■), u d (■)) with P d ( 0) = -Pp Cgs o yielding a larger objective value. There 
exist two other equilibria satisfying SPP, namely the two FCSS 0 , the oligotrophic (-^pcss 0 ’"^FCSS 0 ) 
and the eutrophic FCSS 0 (-^pcss 0 ’ ^FCSS 0 ) 6a). For each of these equilibria there may exist 

a path with P d (0) = Pp CSS o and converging to the oligotrophic/eutrophic equilibrium. To determine 
these paths we consider the corresponding homotopy problems Eq. (25) with x\ replaced by Pp CSS o, 
starting at the equilibrium (Pp£ ss o , Pp£ ss o) and (Ppcss 0 ’ K’css 0 ^ ’ res P ectivel y- 

4.3.1 Comparison with equilibrium solution 

The result of these computations is depicted in Fig. 6. Figure 6a illustrates the “embedding” Eq. (25b). 
The green manifold proceed from the flat eutrophic states -Pp^ggo and the blue manifold from the flat 
oligotrophic states PpQ SS o to the patterned states -Pp Cgs o • During the continuation of k from zero to 
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one the initial states lie in these manifolds. In Figs. 6b and 6c the states and their norm of the solution 
paths for k = 0.5 are shown (green and blue). Additionally the states and norm of the equilibrium 
solution -Pp CSS o are depicted. 

Figures 6d to 6f display the final results for k = 1 also including the figure with the control paths. 
Moreover, Fig. 6g and Fig. 6e reveal that the values of the costates and hence the controls of the three 
solutions with P d { 0) = Pf css o are different. 

Figure 6h show the slice manifolds (see Definition A.l) in the state-costate space. Each marker x 
denotes a continuation step. 

In Fig. 6e the objective values along the slice manifolds are plotted, i.e. the objective values 
are plotted against the (norm) of the initial points. We see that the final solutions for k = 1 both 
yield a higher objective value and the eutrophic solution converging to FCSS °2 dominates the other 
solutions. On the other hand the gradients of the curves suggest that they eventually intersect. Such 
an intersection point characterizes an indifference threshold point (cf. Fig. 3b for the 0D model). 

4.3.2 Detection and Continuation of an indifference threshold point ITP 

To find a possible intersection point of the objective values along the slice manifolds we have to assure 
that the slice manifolds are comparable and have to check if they intersect (see Definition A. 2). The 
slice manifolds depicted in Fig. 6h are, e.g., not comparable, simply because the manifolds 

{■^FCSS 0 — “iH^PCSS 0 ~ ^FCSS°) ' ai ^ 

{■^FCSS 0 (1 — a2 )(-Pp CSS 0 - Pp’ css o) '■ Q 2 e *}■ 

are different. See Fig. 6a), where the green and blue curves depict (parts) of these manifolds. Only 
if slice manifolds are comparable and have a non-empty intersection this means that the solutions 
corresponding to the intersection points of the slice manifold start at the same initial states and hence 
their objective values can be compared. 

Anyhow, since the solution starting at Pf css o and converging to FCSS °2 is known we can start 
the homotopy BVP (25) with 

A ( 0 ) = Pp CSS o + (1 - K )iPpcss° ~ ^PCSS°)- 

In that case the manifolds 

i^FCSS 0 ~ ai )(^PCSS° ~ ^FCSS°) ' ai G 
{-^PCSS 0 + (1 ~ k )(-^fcss° ~ ^PCSS°) : K e ^1- 

trivially coincide and the slice manifolds are comparable, see Fig. 7a. Moreover the slice manifolds 
intersect and and the corresponding objective values curves intersect at an indifference threshold point 
P? A (cf. Fig. 7d). 

In Fig. 7b and Fig. 7e the corresponding state Pf a {-) and control paths u d D (-)) are shown. 

Repeating the same procedure for the three two FCSS 0 and PCSS 0 we find that PCSS 0 is not 
optimal and again we find a further indifference threshold point Pf 2 - The according solutions are 
depicted in Fig. 7c and Fig. 7f. 

It is an interesting question to see how we can find indifference threshold points “between” Pf j 
and Pf 2 • Let us recall that we found the indifference threshold points by the intersection of two n 
dimensional manifolds (with n = N + 1 the number of states P d ) in the n + 1 dimensional space 
(P d 1 J(P d ( 0))). Generically this yields an n — 1 = N dimensional manifold. Already for the dis¬ 
cretization N = 50 its dimension is too large to recover the entire indifference threshold manifold (in 
the original PDE problem the dimension actually increases to infinity). Anyhow, we can e.g. search 
for indifference threshold points along the linear connection Pf 2 + (1 — n)(Pf 1 — Pf 2 ). The result is 
depicted in an animation embedded into Fig. 7e. The according BVP for the numerical solution of 
this problem is presented in Appendix A.2.2. 
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(d) State paths at k = 1 


(e) Control paths at k = 1 


0 20 40 60 60 100 

t 

(f) Norm of the states at ft = 1 



Figure 6: This figure presents the steps of the continuation process to find solutions starting at 
P d ( 0) = Pp CSS o and converging to the oligotrophic and eutrophic FCSS 0 . The colors refer to the 
eutrophic (blue), oligotrophic (green) and patterned (magenta) solutions. In (a) the manifolds of 
initial distributions, that are passed through the continuation, are depicted. For the continuation 
parameter k, = 0.5 the state paths and corresponding norms are shown in (b) and (c). The final 
results are illustrated in (d) (state paths), (e) (control paths), (f) (costate paths), and (g) (norms). In 
(h) the corresponding slice manifolds are shown in the state-costate space. The objective values for 
solutions of the slice manifolds in (i) show that the path converging to the oligotrophic equilibrium is 
optimal among all solutions that start in P d ( 0) = Pp CSS o- 
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(d) Crossing of objective values utioi (f) Control paths of Skiba solution 

Figure 7: This figure shows the detection of a patterned Skiba distribution and its continuation to a 
different Skiba distribution on the Skiba manifold. To receive the animation files associated to the 
panels (a), (b) and (e) please contact the author. 
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4.4 Second scenario 

In the first scenario with b = 0.65 we only found the eutrophic and oligotrophic equilibria being FOSS. 
This is somehow analogous to the 0D model, where the CSS does not appear in the optimal system 
(see Fig. 3). 

The result for the 0D model in the second scenario shows that beside the eutrophic and oligotrophic 
CSS 0 the unstable node (CSS - ) appears as limit of the regions of attractions for the two CSS 0 . Thus, 
we can expect that FCSS - is optimal in model (20). Therefore, also PCSS - cannot be excluded to 
be optimal. 

We therefore analyze two different cases with c = 3.5 where none of the PCSS satisfy SPP and 
c = 3.0825, for which one of the PCSS satisfies SPP. For the original shallow lake model (16) these 
two cases are qualitatively the same (cf. Fig. 4). 

4.4.1 Patterned equilibrium not satisfying SPP 

For c = 3.5 we try to find solutions that start at the states of the FCSS - {Pp CSS ~) and one of the 
PCSS - and converge to the eutrophic (oligotrophic) FCSS. 

In Fig. 8 the main results of this analysis are depicted. The example depicted in the first column 
(a,d and g) is analogous to the case in the 0D model shown in Fig. 4. Thus, it is not possible 
to find a solution starting at the initial states of the FCSS - and converging to the eutrophic or 
oligotrophic FCSS. Instead, during the continuation process the initial states of Pp CSS - are approached 
but cannot be reached. In Fig. 8a the phase portrait of the final result is plotted, which is analogous 
to the state-costate space in 0D (cf. Fig. 4a). In Fig. 8c we see that the objective value for the 
homogeneous initial distributions is continuous (cf. Fig. 4b). Consequently, for spatially homogeneous 
initial distributions, the optimal path is unique, where Pp CSS - separates the regions of attractions 
for the eutrophic (oligotrophic) equilibrium PpQ SS o (Pp° 0 ). The optimal state paths for solutions 
starting exactly at Pp CSS - (the equilibrium solution, black) and in the near vicinity (blue and green) 
are depicted in Fig. 8d. Repeating these steps for each of the PCSS, two examples are depicted in the 
last two columns of Fig. 8, yields that each of these equilibria is optimal, i.e. are POSS. Since none of 
the PCSS satisfy SPP these equilibria and their stable manifolds separate the regions of attractions 
of the FOSS. 

At this point a few questions remain unsettled. 

• Is the defect of an equilibrium not satisfying SPP constant for state discretizations that are fine 
enough? 

• What does this mean for the original PDE problem? 

• Can we say that equilibria not satisfying SPP and their stable manifolds separate the regions of 
attractions for the solutions satisfying SPP? 

• What is the state space of the PDE problem? 

4.4.2 Patterned equilibrium satisfying SPP 

In this section we numerically check, whether the unique patterned equilibrium PCSS 0 for c = 3.0825 
is POSS(cf. Fig. 5b). For that reason we have to show that there exists no other solution path of 
the canonical system Eq. (21) that starts at Pp CSS o yielding a larger objective value. The only other 
candidates are stable paths of the eutrophic and oligotrophic FCSS 0 . The result of the numerical 
comparison for the oligotrohic versus the patterned equilibrium is depicted in Fig. 9d and Fig. 9e. 

To find a feasible path (Pf(-, k), Xf (•, ac)) that satisfies P d ( 0,1) = -Pp C g S o and 

lim (Ff(f, k), A f(t, «)) = (J*£ ss Odessa) with [0,1] 
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we solve the homotopy problem Eq. (25), starting with the constant equilibrium solution (P d '° ss0 , P d £ ss o)- 
The continuation process revealed that it is not possible to find a feasible path for k = 1, instead some 
value kq < 1 was approached. The last computed path (Pf(-, kq), Xf (■, kq)) is shown in Fig. 9a together 
with the corresponding slice manifold (dashed black). 

Next we repeated the procedure for the reversed homotopy problem, starting with the constant 
patterned solution (Pp CSS o, ^pcss°) anc ^ trying to find a feasible path (Pj (■, 1 — k), 1 — ft)) that 

satisfies P d ( 0,0) = PpQ Ss o and 

t lim (P|(t, 1 - ft), \%(t, 1 - ft)) = (Pp CS go, Ap CSS o) with ft G [ 0 , 1 ]. 

Again the continuation process revealed that it was not possible to find a feasible path for k = 0, 
instead some value approximately 1 — fto was approached. This solution (Pf (■, 1 — ft), A^-, 1 — ft)) is 
represented in Fig. 9b by the blue solution path and black slice manifold. 

The two last solution paths from both continuation processes suggest that there exists a separating 
manifold for the regions of attractions of the oligotrophic FCSS 0 and PCSS 0 . A possible candidate 
for this separating manifold is the stable manifold of the PCSS - with defect —1 (see the o in Fig. 9c). 

To test this conjecture we solved the homotopy problem Eq. (29) for defective equilibria. For 
we took 0 , kq) (the initial states of the last continuation step of the first homotopy problem) 
and set V\ := (1,...,1) G l^ 1 , which satisfies the rank condition Eq. (29d). The last solution 
(P 3 (•, 1), A 3 (-, 1)) of this homotopy problem is depicted as dashed blue curve in Fig. 9d and gives a 
strong numerical argument for our conjecture. 

The overall picture, Fig. 9d, suggests that for every e > 0 there exists K\ and ft 2 such that there 
exists solutions (P 1 rf (-, fti), Af(-, fti)) and {P^i', 1 — ft 2 ), Af^-, 1 — k 2 )) of the homotopy problems with 

||(f?(0,/c 1 ),Af(0,K 1 )) - (ff(0,l),Ag(0,l))|| a < £ 

and 

||(-f *2 (0,1 - « 2 ), Ag(0,1 - « 2 )) - (P 3 d (0,1), A 3 ( 0 ,1 ))|| 2 < £ 

or even stronger 

ll(P?('. «i). 4(; «i)) - (Pi (•■ l),Ag(-, l))|| Ia < £ 

and 

W(P d (; 1 - « 2 ), Ag(-, 1 - « 2 )) - (P 3 d (-, 1), Af (•, l))|U a < e. 

Plotting the objective values evaluated along the solutions of the corresponding slice manifolds shows 
that the objective function is continuous in the vicinity of P d (0,1), see Fig. 9e. An analogous result 
holds for the comparison of the stable paths converging to the eutrophic FCSS 0 and the PCSS 0 . This 
proves the optimality of FCSS 0 and PCSS 0 as well. Also according to the case c = 3.5 the stable 
manifolds of defective equilibria separate the regions of attractions of the equilibria satisfying SPP. 

Figure 9f shows (part of) a phase portrait in the state-costate space for c = 3.0825. The sub¬ 
scripts of the equilibria denote the defect of the according equilibrium. Thus, there exist two FCSS 0 
(PpQ Sg0 and Pp°, ggo), and one PCSS 0 {Pp CSS o)- Additionally a few paths are plotted that converge to 

^fcss 0 ’ ^FCSS 0 ’ an< ^ ^PCSS 0 ( so ^d blue) and to PCSS - with defect —1 (dashed blue). The specific case 
for solutions that converge to the oligotrophic equilibrium (PpQ Ss oi -^fcss°) an< ^ Patterned equilibrium 
(.Ppcgso 5 -^pcss 0 ) separately illustrated in Fig. 9d and Fig. 9e. 

5 Conclusion 

The purpose of this article is the presentation of a numerical framework, that allows us to numerically 
experiment with ID spatially distributed optimal control problems. Numerical experiment is meant 
in the sense of Heaviside, who claimed mathematics as an experimental science, cf. feaviside [1893]. 
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(g) Objective value along slice manifolds (h) Objective value along slice mani- (i) Objective value along slice manifolds 

folds 


Figure 8: This figure depicts the numerical proof for the optimality of FCSS - and PCSS - at c = 3.5, 
i.e., the optimality of the equilibria that do not satisfy SPP, exemplified for three equilibria. In the 
first row (a)-(c) the slice manifolds for the according continuation processes are depicted in the norrned 
state-costate space. The second row (d)-(f) shows the state paths for the solutions starting at and 
near the FCSS - and PCSS - , respectively. The last row (g)-(i) illustrates that the objective function 
is continuous in the vicinity of the constant equilibria solutions for the FCSS - and PCSS - . Therefore 
the equilibria not satisfying SPP are optimal as well. 
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(e) Objective value (f) Phase portrait 


Figure 9: In (f) some of the solutions paths and equilibria in the state-costate space are depicted. 
The solutions shown in (d) and (e) are the result of the continuation processes, when we tried to find 
a solution starting at the states of POSS and converging to the oligotrophic FOSS, and vice versa. 
The continuation process approached a state lying on the stable manifold (dashed blue) of the PCSS - 
with defect —1. Thus, for initial states coinciding with the states of the stable manifold with defect 
— 1 it is optimal to converge along the stable manifold to the POSS with defect —1. In the vicinity 
of these initial states it is optimal to converge either to the oligotrophic FOSS or the POSS satisfying 
SPP. To receive the animation file associated to the panel (d) please contact the author. 
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We were thus able to find numerical evidence for the occurrence of (indifference) threshold distri¬ 
butions in distributed models. Moreover, this points towards a possible generalization of the saddle 
point property, which in case of multiple canonical steady states (CSS) is intimately connected to 
the existence of multiple optimal solutions. Even though this simple FDM could also be applied to 
spatially 2D models, such an approach would immediately become numerically intractable. Therefore, 
it is an intermediate step to the method of a finite element discretization, that is presented in Irass 
and Ueckcr [201 ] and extended in [2015]. 

Different directions for further research result from the presented approach. The obvious next step 
is the previously mentioned application of FEM discretization. This is realized as an add-on toolbox 
(p2pOC) to the MATLAB package pde2path, which is a numerical tool for continuation and bifurcation 
in 2D elliptic systems, cf. Uecker et al. [2014]. 3 

A main drawback of the actual approach is difficulty to handle (inequality) constraints. For a 
correct usage of the used BVPsolver the rhs of the dynamics has to at least continuously differentiable. 
This property is violated, when constraints become active. We solved that problem by considering 
different arcs of the solution path, where each of the arcs satisfied the differentiability condition. For 
non-distributed models this is in general a practicable way but hardly to realize for distributed models. 
Simply because for each transition from one arc to the other the according switching conditions have 
to be stated. This ansatz quickly becomes intractable for the high dimensional discretized system. 

Therefore, we will put effort in the development of a solver based on finite element discretization 
and conjugate gradient method combined with a continuation step. 


A Numerical method implemented in OCMat 

In this section we formulate the basic method for the calculation of paths that converge to an equi¬ 
librium of the canonical system. 


A.l The core problem 

The core problem that has to be solved is the following. Given an equilibrium (x, A) and an initial 
distribution xo we want to find a path (x(-),A(-)) satisfying Eq. (3) together with the boundary 
conditions 

x(0) = xo and lim (x(t), A(t)) = (x, A). ( 22 ) 

£—>■00 

The dimension of the according eigenspace 

0 < dim E s (x, A) = n s < n. 


Then, the boundary condition at infinity in Eq. (22) can be approximated by the so called asymptotic 
boundary condition [cf. .entini, 1978, Lentini and Keller, 1980] 


P T 




0GR 2n-n % ^ GR 2nx(2n-n s ) &nd n±E »( £> X) 


(23) 


with T > 0 large enough. 

For a compact notation we introduce 


X : = 



and Eq. (21) is written as X(t.) = F(X(t)). 


3 It can be downloaded for free from http://www.staff.uni-oldenburg.de/hannes.uecker/pde2path. 
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Then the previous BVP writes as 


(24a) 

(24b) 

(24c) 


X(t) = TF(X(t)), t G [0,1] 

I 2 T (V - X(l)) = 0 . 

If X satisfies SPP then n s = n and Eq. (24b) simplifies to X^ 1 , (0) = xq. To keep notation simple 
we assume that the coordinates of x are sorted, such that (*i,... ,i n „) = (1, ■ ■ • ,n s ) =: (n s ). Then 
Eq. (24b) can be rewritten as 


X M {0)=x ( o ns) (24b’) 

In general BVP (24) cannot be solved analytically, hence numerical methods have to be applied. These 
numerical methods request some initial function X(-). Such an initial function need not satisfy the 
BVP (24) but, depending on the problems properties, it has to a more or less good approximation. 
What can we do if such an initial function is not at hand? 

A.2 Embedding into a homotopy problem 

Given that a solution Y(-) of Ecp (24) with Y^ n \ 0) = =4 is available we can embed Eq. (24) 

into the according homotopy problem 

X(t) = TF(X(t),n), t G [0,1] (25a) 

I w (0) = x ( 0 n) + (1 - K)(xj n) - 4 n) ) (25b) 

52 T (A - A(l)) = 0. (25c) 

Then Y (•) solves Eq. (25) for k = 0 and k = 1 yields a solution of BVP (24). OCMat solves the homotopy 
BVP (25) using arclength continuation [Kuznetsov, 1998, Allgower and Georg, 2003]. Thus, at each 
homotopy step i > 0 with previous solution (A(j_ 1 )(-), «j_i) the BVP (25) is solved for (Vqq(-),^i) 
together with the additional equation 


[ ~ x (i-i)(t)dt + (m - = 0 (25d) 

Jo 

where (Vfj_ij(-), V(,-_i) )K ) satisfies the linearized BVP 

V(i-i)(t) = TEx(V( i _ 1 ) (t), /r)V (i _ 1 ) (t), t G [0,1] (25e) 

V ( S)(0) = V {i _ 1)tK (x^ - 4 n) ) (25f) 

H T V( 1 ) = 0. (25g) 


The solution (Vo_i)(-), Vu_i^ K ) of BVP (25e)-(25g) is called the tangent solution at step i — 1. In 
the actual OCMat implementation the BVP (25a)-(25g) is discretized providing different discretization 
schemes, relying on the two native MATLAB BVP solvers bvp4c, bvp5c [cf. vierzenka and Shampine, 
2001, 2008], and the adapted solver bvp6c [cf. Hale, 2006, Hale and Moore, 2008]. The discretized 
tangent is computed at each Newton step. This is a specific arclength continuation, called Moore- 
Penrose continuation, cf. 4uznetsov [1998]. 

Subsequently we introduce some terminology to make a clear distinction between the stable paths 
and the set of initial points computed during a continuation process. 
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Definition A.l (Slice Manifold). Let X(-,k(s)), s 6 /Cl with I a non-empty interval and k(-) G 
C°(I, M) be a solution of Eq. (25) for every s G I. Then 

S{X,x^\xf s \n{-)) := {X(0,«(a)) : s € 1} (26) 

is called the slice manifold along for X and k(-). 

Definition A.2 (Comparable Slice Manifolds). Let X satisfy SPP and S(Xj, Xq, Kj(sj)), Sj G 
Ij , j = 1,2 be two slice manifolds along Xq,x{ for Xj and Kj(-) with Xq x\. j = 1,2. Then 
S(Xj, Xq, x \, Kj{sj)). j = 1,2 are called comparable iff 

{xj + (1 - oti)(x\ - Xg) : Qi G 1} = {xq + (1 - a 2 )(xj - Xq) : a 2 el} (27) 

holds. 

If comparable slice manifolds satisfy 

{xq + (1 - Ki(s))(x\ - Xq) : s G h} n {xg + (1 - re 2 (s))(x? - Xq) : s G I 2 } ^ 0 (28) 

it is said that the slice manifolds are intersecting. 

Remark A.l. A slice manifold is a linear cut through the stable manifold. At the intersection of 
two different slice manifolds the cuts are given for the same (initial) states xo- Hence the according 
paths are (different) solution candidates for the optimal control problem with x(0) = xq. For one state 
autonomous optimal control problems with the stable path (x(-), A(-)) converging to (x, A) 

5(x 0 ,x(r),x,id [0il] ) = {(x(f),A(f)) : t G [0,T]}. 

Thus, the orbit of the stable path coincides with the stable manifold. Moreover two slice manifolds for 
different saddles are trivially comparable. 

There are good reasons to consider BVP (25) instead of BVP (24). For an arbitrary initial point 
xo it is often hard to provide a “good” guess of an initial function for BVP (24). Since in general the 
solution of an BVP is not unique it may not be guaranteed that a computed solution is the searched 
for solution. 

On the other hand the equilibrium solution trivially satisfies BVP (24) for the initial point x. Hence 
the homotopy BVP (25) can be started with an exact solution. Then the existence of a unique solution 
is guaranteed by the implicit function theorem as long as some rank condition is satisfied. A careful 
inspection of the linearization of BVP (25), which is a byproduct of the arclength continuation, then 
yields important information about the behavior of the solution paths. Moreover, these intermediate 
solution paths can be used to find, e.g. an indifference threshold point. 

A.2.1 The stable manifold for an equilibrium not satisfying SPP 

In the previous section we assumed that X satisfies SPP, i.e., dimE s (J) = n. Next we adapt BVP (25) 
for the case dimE s (J) = n s < n. Therefore, we assume that X is hyperbolic with dimE“(J) = n u 
and hence 2 n = n s + n u . Furthermore, two points x) ' in the state space are fixed and n — n s vectors 
Vi G are chosen. Then, the according BVP for the calculation of stable paths converging to X 
becomes 

X(t) = TF(X(t),n ), t G [0,1] (29a) 

n—n s 

-V (ti) ( 0) = X^ n) + (1 - Ko){Xo l] - X^ n) ) + ^2 K i y i ( 29b ) 

2=1 

n T (X - X(l)) = 0 G M n “ (29c) 
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with 

rank ((xj, n) - x^) v x ••• v n - n )j = n - n s + 1 (29d) 

OLE S {J) = 0 and lie M 2nxn + 

Let us assume that a solution Y(-) for Kj = 0, j = 0,..., n — n s is known. Then Eq. (29b) can be 
interpreted as the condition that we search for a solution along the direction . Additionally 

we have to take care letting enough freedom, guaranteed by 29d, to start at the stable manifold. 

A.2.2 Continuation of an indifference threshold point 

A useful relation between the Hamiltonian and objective value is given for a model (1) with p > 0 
and finite objective value [cf. Michel, 1982]. Then we find for any solution x(-),\(-) of the canonical 
system Eq. (3) 

J(so) = -ft(*(0),A(0)) (30) 

P 

where Ti is defined according to Eq. (3d) and the bar is omitted. 

Let Yi, i = 1, 2 be two CSS of Eq. (3) and x} and xj be two distinct indifference threshold points 
of model (1). Furthermore, let -^ 1 , 2 (•) be two solutions corresponding to x} and ^ 3 , 4 (•) two solutions 
corresponding to x’j. To continue the indifference threshold point from xj to x'j we solve the following 
homotopy problem 


X 1 (t) = T 1 F(X 1 (t)), t e [ 0 , 1 ] 

(31a) 

X 2 (t) = T 2 F(X2(t)), t € [0,1] 

(31b) 

A{ n) (0) = Xi; n) (0) <E M n 

(31c) 

H(Ai( 0 ))-H(X 2 ( 0 )) = 0 eK 

(31d) 

H 7 (Ti-Ai(l)) = 0 GM n 

(31e) 

^(y 2 -i 2 (i)) = oer 

(31f) 

Aj n) (0) = xj + (1 - «i)(x} - xj) + k 2 V g M n 

(31g) 

a\V + 02 (xj — xj) = 0 and 04 + 02/0 


Lb+E s (Jj), * = 1,2. 



Equations (31a) and (31b) denote the dynamics for the two distinct paths, starting at the same initial 
states, Eq. (31c). The truncation times T\ > 0 and To > 0 may be chosen differently. The two paths 
Xl(-) and X 2 (-) yield the same objective value which according to Eq. (30) can be stated as Eq. (31d). 
Equations (31e) and (31f) denote the asymptotic boundary conditions for the path Ai(-) converging 
to Y\ and A^-) converging to Y 2 , respectively. Finally Eq. (31g) specifies the continuation in the state 
space, where the first part xj + (1 — ki)(x) — xj) describes the change into the direction to the target 
xj and K 2 V is a correction term, since the stable manifold has one dimension less than the state space. 

Counting the number of unknowns and equations, we find 4n + 2 unknowns, two times the states 
and costates A,;(•) and the two free parameters Kj, i = 1,2 and 4n + 1 equations. Moreover Eq. (31) 
is solved for 0,0) and (^ 3 , 4 (-), 1, 0). Thus we can start a continuation process starting with 

one of these previously detected solutions. 
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B The usage of OCMat 


In this section the basic steps for the numerical analysis of model (20) are explained in detail. This 
enables the user to reproduce the presented results and learn the basic commands and structure of 

OCMat. 

B.l The initialization file 

To get results that are comparable smooth in the spatial dimension as the discretization used in Grass 
and Ueckcr [20! ] we choose IV = 51. In the syntax logic of OCMat we have to provide an initialization 
file consisting of IV + 1 ODEs, the entry of IV + 1 state and control variables Pi, Ui, i = 0,..., N and 
the appropriate objective function. Doing that by hand can become a boring and error-prone task. 
We therefore wrote a MATLAB file ’makeinitf ile ’ that generates the initialization file, providing 
IV 4 

function makeinitfile(N,modelname) 
if nargin==l 

modelnaine=’shallowlakeline’; 

end 

dotPxi=’uxi-b*Pxi+Pxi~2/(1+Pxi~2)+D*(N/2/L)~2*(Pximl-2*Pxi+Pxipl)’; 
intvar =’’; 
for ii=0:N 

ode{ii+l}-=[’DPx’ num2str(ii) ’ = ’ strrep(strrep(strrep(strrep(strrep(dotPxi,’i’,num2str( 
ii)),[’x’ num2str(ii) ’ml’],!’*’ num2str(ii-1)]),Ex’ num2str(ii) ’pi’],Ex’ num2str(ii 
+1)]) ,’Px-1’,’Pxl 1 ), [’Px’ num2str(N+l)],[’Px’ num2str(N—1)])]; 
ode{ii+ll=char(simple(sym(ode{ii+l}))) ; 
if ii>0 

controlvar=[controlvar Dux’ num2str(ii)]; 
statevar=[statevar ’,Px’ num2str(ii)]; 
if ii<N 

intvar=[intvar ’+log(ux’ num2str(ii) ’)-c*Px’ num2str(ii) ’~2’ ]; 

else 

intvar=[intvar ’+(log(uxO)-c*PxO~2+log(ux’ num2str(ii) ’)-c*Px’ num2str(ii) ’ 

~ 2 )/ 2 ’ ]; 
end 

else 

controlvar=[ , ux’ num2str(ii)]; 
statevar=[’Px ) num2str(ii)]; 

end 

end 

intvar (1) = [] ; 

intvar=[’(’ char(collect(simple(sym(intvar)),’c’)) 

par={ ’ rho : :2*pi/0.44 ) , [’N: : ’ num2str(N)] >; 

initfilefid=fopen([modelname ’ .ocm’],’w’); 

fprintf (initf ilef id, ; Type\nstandardmodel\n\nVariable\nstate : :7,s\ncontrol: :7„s\n\ 
nStatedynamicsXn',statevar,controlvar); 
for ii=l:length(ode) 

fprintf (initf ilef id, ’ode: :7 0 s\n’ ,ode{ii}) ; 

end 

fprintf(initfilefid, 1 \nObjective\nexpdisc::rho\nint: ^/osXnXnParameterXn 1 ,intvar); 
for ii=l: length, (par) 

4 If this FDM approach turns out to be an important tool by itself, this step could be directly implemented within 
OCMat. 
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fprintf (initf ilef id, ’"/.sW’ ,par{ii}) ; 

end 

fclose(initfilefid); 

As an example we call 
>> makeinitfile(5,’shallowlakelinetest’) 
at the MATLAB workspace and the file ’shallowlakelinetest.ocm’ is generated 
Type 

standardmodel 

Variable 

state::PxO,Pxl,Px2,Px3,Px4,Px5 
control::uxO,uxl,ux2,ux3,ux4,ux5 

Statedynamics 

ode::DPxO=uxO-b*PxO+Px(T2/(1+Px0'2)+D*(N/2/L)~2*(Pxl-2*PxO+Pxl) 
ode::DPxl=uxl-b*Pxl+Pxl"2/(l+Pxl~2)+D*(N/2/L)~2*(Px0-2*Pxl+Px2) 
ode::DPx2=ux2-b*Px2+Px2~2/(1+Px2~2)+D*(N/2/L)~2*(Pxl-2*Px2+Px3) 
ode::DPx3=ux3-b*Px3+Px3~2/(1+Px3~2)+D*(N/2/L)~2*(Px2-2*Px3+Px4) 
ode::DPx4=ux4-b*Px4+Px4~2/(1+Px4~2)+D*(N/2/L)~2*(Px3-2*Px4+Px5) 
ode::DPx5=ux5-b*Px5+Px5~2/(1+Px5'2)+D*(N/2/L)~2*(Px4-2*Px5+Px4) 

Objective 
expdisc::rho 

int::((-Pxl~2-Px2~2-Px3~2-Px4~2-l/2*Px0~2-l/2*Px5~2)*c+log(uxl)+log(ux2)+log(ux3)+log(ux4) 
+l/2*log(ux0)+l/2*log(ux5)) 

Parameter 
rho::0.03 
b::0.65 
c : : 0.5 
D : : 0.5 

L::2*pi/0.44 
N: :5 

This file is placed at the actual MATLAB directory, to make it visible for OCMat it has to be moved 
to the folder ocmat\model\initf iles. After that the initialization process of OCMat can be started 

>> ocStruct=processinitfile (’ shallowlakelinetest’); 

ocmat\model\usermodel\shallowlakelinetest does not exist. Create it? (y)/n: y 
ocmat\model\usermodel\shallowlakelinetest\data does not exist. Create it? (y)/n: y 
ocmat\model\usermodel\shallowlakelinetest\data is not on MATLAB path. Add it? (y)/n: y 
ocmat\model\usermodel\shallowlakelinetest is not on MATLAB path. Add it? (y)/n: y 
>> modelfiles=makefile4ocmat(ocStruct); 

>> moveocmatfiles(ocStruct,modelfiles) 

For the following analysis we use the initialization file shallowlakelinecoarse with N = 51. Thus, 
the previous initialization commands have to be repeated with shallowlakelinecoarse instead of 
shallowlakelinetest. Note that even though the discretization parameter N appears as a parameter 
of the model shallowlakelinecoarse it must not be changed. The initialization steps can be very 
time consuming, depending on the chosen value of N and the computer capacity. For A^ = 51 the 
initialization, on a PC with the specifications Intel Core i7-3820 CPU@3.60GHz, Windows 7, MATLAB 
7.6.0, took around half an hour. 

Subsequently we will also make use of the 0D shallow lake model. Therefore, we have to run 
through the initialization process for the model shallowlake 0 as well. 

5 The according initialization file shallowlake. ocm can be found in the folder ocmat/mode/initf iles. 
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>> ocStruct=processinitfile(’shallowlake’); 

ocmat\model\usermodel\shallowlake does not exist. Create it? (y)/n: y 
ocmat\model\usermodel\shallowlake\data does not exist. Create it? (y)/n: y 
ocmat\model\usermodel\shallowlake\data is not on MATLAB path. Add it? (y)/n: y 
ocmat\model\usermodel\shallowlake is not on MATLAB path. Add it? (y)/n: y 
>> modelfiles=makefile4ocmat(ocStruct); 

>> moveocmatfiles(ocStruct,modelfiles) 

B.2 Detection and continuation of equilibria 

For the calculation of the FCSS we use the equilibria of the shallowlake model. A detailed explanation 
of the commands can be found in the OCMat manual. 

>> mO=stdocmodel( 1 shallowlake’); 

>> mO=changeparametervalue(mO,’b,c’,[0.65 0.5]); 

>> ocEP0=calcep(m0);b=isadmissible(ocEPO,mO,[],’UserAdmissible’);ocEP0(~b) = [] ; 

The values of the equilibria are used to generate the values of FCSS. 

>> m=stdocmodel(’shallowlakelinecoarse’); 

>> N=parametervalue(m,’N’); 

>> y0=[ocEP0fl}.y([1 end]) ocEP0{2}.y( [1 end]) ocEP0{3},y( [1 end])]; 

>> Y= [y0( [ones(N+l,1);2*ones(N+l,1)],1) y0([ones(N+l,1);2*ones(N+l,1)],2) y0( [ones(N+l,1) 

;2*ones(N+l,1)],3)] ; 

» Y( [N+2 2*N+2],:)=Y([N+2 2*N+2],:)/2; 

>> opt=setocoptions(’EQ ’,’TolFun’,le-12,’MaxFunEvals’,50000, 1 Maxlter’,50000); 

>> ocEP=calcep(m,Y,[],opt);b=isadmissible(ocEP,m);ocEP(~b)=[]; 

>> [b dfct]=isspp(m,ocEP{:}) 
b = 

10 1 

df ct = 

0-5 0 

Thus, the first and third equilibrium satisfy SPP, the second equilibrium has defect —5. For the 
subsequent step an adapted version of CL_MATC0NT has been used. 

>> opt0=setocoptions(’MATC0NT’,’MaxNumPoints’,750,’MaxStepsize’,le-1,’Backward’,0,’ 
IgnoreSingularity’,2,’0CC0NTARG’,’CheckAdmissibility’,’off’); 

>> optl=setocoptions(optO,’MATC0NT’,’Backward’,1); 

>> contpar=’b’;epidx=l;[xO vO sO fO hO]=contep(m,ocEP{epidx]-,contpar,optO); 
first point found 

tangent vector to first point found 
elapsed time = 107.2 secs 
npoints curve = 750 
>> store(m,’modelequilibrium’) 

During the continuation of the first equilibrium four branching points are detected. These branching 
points are the initial solutions for heterogeneous equilibria. Repeating the previous steps for the third 
equilibrium epidx=3 in both continuation directions yields the eutrophic arc. On this arc no branching 
point exists. With store(m, ’modelequilibrium’) the results of the continuation process are stored 
into the results of model m. 

>> contpar=’b’;epidx=3; [xO vO sO fO hO]=contep(m,ocEP{epidx},contpar,optO); 
first point found 

tangent vector to first point found 
elapsed time =72.9 secs 
npoints curve = 750 
>> store(m,’modelequilibrium’) 

>> contpar=’b’;epidx=3;[xl vl si fl hl]=contep(m,ocEP{epidx},contpar,optl); 
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first point found 

tangent vector to first point found 
elapsed time =76.2 secs 
npoints curve = 750 
>> store(m,’modelequilibrium’) 

The calculation of the arcs emanating from the branching points are exemplified for the second branch¬ 
ing point. 

>> n=3;ml=changeparametervalue(m,’b’,[x0(end,s0(n).index)]); 

>> ocEPl=calcep(ml,x0(l:end-1,s0(n).index),0,opt); 

>> optO=setocoptions(optO,’MATCONT’,’MaxNumPoints’,1500,’MaxStepsize’,2e-l,’CheckClosed’ 
,50); 

>> epidx=l;[xbpO vbpO sbpO fbpO]=contbp(ml,ocEPl{epidx},sO(n),0.01,opt0); 
first point found 

tangent vector to first point found 

label = BP, x = ( 0.619422 0.619661 0.620389 ... -5.917192 -5.917341 -2.958687 0.721718 ) 

label = BP, x = ( 0.554176 0.554239 0.554455 ... -6.412585 -6.407783 -3.203096 0.706011 ) 

label = BP, x = ( 0.398031 0.398863 0.401801 ... -9.309244 -9.141784 -4.542734 0.626785 ) 

elapsed time = 249.5 secs 

npoints curve = 1500 

>> store(m,’modelequilibrium’) 

During the continuation starting from the second branching point, i.e. the bifurcation point stored in 
the structure sO at index n=3, further branching points are detected. The branches emanating from 
these PCSS can be computed in an analogous way. Repeating the previous steps for n=2,4,5 finally 
yields all branches of FCSS and PCSS, see Fig. 5a. 

Using the solutions of the previous calculations can now be used to find the all FCSS and PCSS 
equilibria for a specific value of b, e.g. b = 0.65 

>> matRes=matcontresult(m);counter=0;b=0.65; 

>> for ii=l:length(matRes), ... 

x0=[matResfii}.ContinuationSolution.y;matRes{ii}-.ContinuationSolution.userinfo. 
varyparametervalue]; ... 

idx0=cont2idx(x0(end,:),b); ... 
for jj=l:length(idxO); ... 

counter=counter+l;ocEP=calcep(m,xO(l:end-1,idx0(jj)),0,opt); ... 
store(m,ocEP);end,end 

For b = 0.65 there exist 17 equilibria. 6 

We will give a two examples for the calculation of a stable path to an equilibrium that satisfies 
SPP and that does not satisfy SPP. 

B.3 Saddle path calculation 

One of the main tasks of OCMat is the calculation of a stable path converging to an equilibrium of 
saddle-type. The computation is done by solving the homotopy problem Eq. (25) or Eq. (29). We 
start with a problem of the first type. 

B.3.1 Stable path when SPP is satisfied 

As an example we compute, for the parameter values b = 0.65 and c = 0.5, the stable path that starts 
at the states of a defective PCSS (ocEP{9>) and converge to the eutrophic FCSS (ocEP{17}). 7 First 

6 With load(m, ’"/,1.2f ’ , [] , ’b,c’) the model where these equilibria and other numerical results are already stored 
can be loaded into the MATLAB workspace. 

7 The order refers to the results stored in the file shallowlakelinecoarse_b_0.65_c_0.50.mat. 
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we load the model with the stored data of the equilibria. 

>> m=stdocmodel('shallowlakelinecoarse’);m=changeparametervalue(m,’b,c’,[0.65 0.5]); 

>> load(m, '°/„1.2f ' ,1, 'b,c') ; ocEP=equilibrium(m) ; 

Next we determine the flat equilibria satisfying SPP and their size of the states. 

>> idx=find(isflat(m,ocEP{:}) & isspp(m,ocEP{:})) 
idx = 

1 17 

>> P=state (m, ocEP-fl}) ;P(1) 
ans = 

0.4530 

>> P=state(m,ocEP{17});P(1) 
ans = 

1.4370 

To start the continuation process we change some of the default options and call the initialization 
function initocmat_AE_EP. The truncation time T is determined by 


inin ?eE s|Re£| ’ 


where To usually is set to 10. 

>> opt=setocoptions(’0CC0NTARG','MaxStepWidth',1,'InitStepWidth’,5e-l, 1 CbeckAdmissibility 1 , 
'off 1 ,'SBVP0C','MeshAdaptAbsTol’,le-4,'MeshAdaptRelTol',le-3,'GENERAL',’ 
TrivialArcMeshNum',20); 

>> eval=real(eig(ocEP{17}));eval(eval>0)=[];T=10/min(abs(eval)); 

>> sol=initocmat_AE_EP(m,ocEP{17},1:N+1,ocEP{9}.y(1:N+l),opt,'TruncationTime’,T); 

After the initialization process the continuation process is started calling bvpcont. 

>> c=bvpcont(’extremal2ep',sol,[],opt); 
first solution found 

tangent vector to first solution found 

Continuation step No.: 1 
stepwidth: 0.5 
Newton Iterations: 1 
Mesh size: 21 

Continuation parameter: 0.0952864 

Continuation step No.: 15 
stepwidth: 1 
Newton Iterations: 1 
Mesh size: 27 

Continuation parameter: 1.01767 

Target value hit. 
label=HTV 

Continuation parameter=l 

elapsed time = 38.4 secs 
>> store(m,'extremal2ep'); 

>> save(m) 

Finally the result of the continuation is stored in the model-object m, for further use. The save 
command allows to store the entire object, i.e. with the stored results, as a MATLAB data file. 
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B.3.2 Stable path when SPP is not satisfied 

The next example was used for the computation of the second and third homotopy process in Sec¬ 
tion 4.4.2, cf. Fig. 9c. The parameter values are b = 0.55 and c = 3.0825. Thus, in a first step we try 
to find a solution that starts at the states of the flat oligotrophic equilibrium (ocEPfl}) and converge 
to the heterogeneous equilibrium satisfying SPP (ocEP{7}). For that reason we repeat the steps of 
the previous Appendix B.3.1. 

>> opt=setocoptions(opt,’0CC0NTARG’,’MaxStepWidth’,lei, 1 InitStepWidth’ ,5e-l,’ 

MaxContinuationSteps’,30,’SBVPOC’,’MeshAdaptAbsTol’,le-3,’MeshAdaptRelTol’,le-2); 

>> m=changeparametervalue(m,’b,c’,[0.55 3.0825]); 

» load(m,’*/„1.4f’,[],’b,c’) 

The results in the actual model ’m’ are overwritten. Proceed? y/(n) : y 
>> ocEP=equilibrium(m); 

>> eval=real(eig(ocEP{7}));eval(eval>0)=[];T=10/min(abs(eval)) 

>> sol=initocmat_AE_EP(m,ocEP{7},1:N+1,ocEPfl}.y(1:N+1),opt,’TruncationTime 1 ,T); 

>> c=bvpcont(’extremal2ep’,sol,[],opt); 

>> store(m,’extremal2ep’);ocEx=extremalsolution(m);n=length(ocEx); 

During the continuation process we encounter that the detected solution does not end “near” the 
equilibrium, i.e. the used truncation time becomes too short. The reason becomes obvious when 
we have a look on Fig. 9d. The computed solution path approaches a stable path of the defective 
PCSS - . Therefore, the time it takes to stay in the vicinity of this equilibrium increases. To overcome 
this problem we extend the homotopy problem Eq. (25) by letting the truncation time T be a free 
parameter value and adding a further constraint, that guarantees that the solution X() not only ends 
at the the (linearized) stable manifold, but also satisfies 

\\X(l)-X\\ =e, (32) 

with some fixed e > 0. To start this extended continuation process we use the solution after 30 
continuation steps 8 and use the initialization argument ’movinghorizon’ 

>> sol=initocmat_AE_AE(m,ocExfn},[1:N+1],ocEP{l}.y(1:N+l),opt,’movinghorizon’,1) 

>> opt=setocoptions(opt,’0CC0NTARG’,’MaxStepWidth’,le2,’InitStepWidth’,5e0,’ 
MaxContinuationSteps’,70); 

>> c=bvpcont(’extremal2ep’,sol,[],opt); 

>> store(m,’extremal2ep’);ocEx=extremalsolution(m);n=length(ocEx); 

Next we compute a stable path that converges to the defective equilibrium PCSS - . The stable manifold 
of the defective equilibrium (it has defect —1) is of dimension (N + 1) — 1. Taking the initial states 
P d ( 0) of the solution of the last continuation process, in OCMat notation ocEx{n}.y(l:Nl,l)+, there 
exists such that Po = P d (0) + n\v\ with v\ = (1,..., 1) T E M Ar+1 is lying in the TV-dimensional 
stable manifold. Thus, we solve the according homotopy problem. 

>> opt=setocoptions(opt,’0CC0NTARG’,’MaxStepWidth’,lei,’InitStepWidth’,le-1,’SBVP0C’,’ 
MeshAdaptAbsTol’,le-4,’MeshAdaptRelTol’,le-3); 

>> Vl=ones(N+l,1);eval=real(eig(ocEP{8}));eval(eval>0)=[];T=10/min(abs(eval)) 

>> sol=initocmat_AE_EP(m,ocEP{8},[1:N+l],ocExfn}.y(1:N+l,1),opt,’TruncationTime’,T,’ 
freevector’,V1); 

>> c=bvpcont(’extremal2ep’,sol,[],opt); 

The result of the computations can be plotted using OCMat plotting commands. 

>> elf;xcoord=l;ycoord=2;xvar=’spatialnorm’;yvar=’spatialnorm’; 

>> plotcont(m,xvar,xcoord,yvar,ycoord,’contfield’,’ExtremalSolution’,’Index’,[2 4 5]);hold 
on, 

8 We therefore set the option ’MaxContinuationSteps’ ,30. 
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>> plotlimitset(m,xvar,xcoord,yvar,ycoord,’Index’,[1 7 8],’Marker’MarkerSize’,10,’ 
showspp’,1,’showflat’,1); 

>> hold off;figure(gcf) 


B.3.3 Indifference threshold point and manifold 

This example presents in detail the computation of the results from Section 4.3.2. First we load the 
models data, retrieve the equilibria and remove the already stored continuation results. 

>> m=stdocmodel(’shallowlakelinecoarse’);m=changeparametervalue(m,’b,c’,[0.65 0.5] ) ; 

>> load(m, ’7.1.2f ’ , 1, ’b, c ’ ) ; ocEP=equilibrium(m) ; 

>> removeresult(m,’Continuation’); 

The solutions that start at the states of the seventh equilibrium, a PCSS 0 , and converge to the 
eutrophic and oligotrophic equilibrium are computed. 

>> idx=f ind(~isf lat (m,ocEP-[: }) & isspp(m,ocEP-[:})) 
idx = 

4 7 

>> opt=setocoptions(’0CC0NTARG’,’MaxStepWidth’,lei,’InitStepWidth’,5e-l,’CheckAdmissibility 
’,’off’,’SBVPOC’,’MeshAdaptAbsTol’,le-4,’MeshAdaptRelTol’,le-3,’GENERAL’,’ 
TrivialArcMeshNum’,20); 

>> eval=real(eig(ocEP{l}));eval(eval>0)=[];T=10/min(abs(eval)); 

>> sol=initocmat_AE_EP(m,ocEPfl},1:N+1,ocEP{7}.y(1:N+1),opt,’TruncationTime’,T); 

>> c=bvpcont(’extremal2ep’,sol,[],opt); 

>> store(m,’extremal2ep’);ocEx=extremalsolution(m);n=length(ocEx); 

>> eval=real(eig(ocEP{17}));eval(eval>0)=[];T=10/min(abs(eval)); 

>> sol=initocmat_AE_EP(m,ocEP{17},1:N+1,ocEP{7}.y(1:N+l),opt,’TruncationTime’,T); 

>> c=bvpcont(’extremal2ep’,sol,[],opt); 

>> store(m,’extremal2ep’);ocEx=extremalsolution(m);n=length(ocEx); 

Plotting the objective value (Hamiltonian) shows that the eutrophic solution is the optimal solution. 
Next the last eutrophic solution is continued (for ten steps) in direction of the oligotrophic equilibrium. 

>> elf;xcoord=l;ycoord=l;xvar=’spatialnorm’;yvar=’hamiltonian’; 

>> plotcont(m,xvar,xcoord,yvar,ycoord,’contfield’,’SliceManifold’,’Index’,[1 2]);figure(gcf 

) 

>> opt=setocoptions(opt,’0CC0NTARG’,’InitStepWidth’,2.5e0,’MaxContinuationSteps’,10); 

>> sol=initocmat_AE_AE(m,ocEx{2},1:N+1,ocEPfl}.y(1:N+1),opt); 

>> c=bvpcont(’extremal2ep’,sol,[],opt); 

>> store(m,’extremal2ep’);ocEx=extremalsolution(m);n=length(ocEx); 

Comparing the objective value of the solutions for the first and third continuation process yields 
an (heterogeneous) indifference threshold point. Finally the solutions starting at the indifference 
threshold point and converging to the FOSS are computed. 

>> ipt0=findindifferencepoint(m,1,3); 

>> opt=setocoptions(opt,’0CC0NTARG’,’InitStepWidth’,lei,’MaxContinuationSteps’,inf); 

>> eval=real(eig(ocEP{l}));eval(eval>0)=[];T=10/min(abs(eval)); 

>> sol=initocmat_AE_EP(m,ocEPfl},1:N+1,iptO,opt,’TruncationTime’,T); 

>> c=bvpcont(’extremal2ep’,sol,[],opt); 

>> store(m,’extremal2ep’);ocEx=extremalsolution(m);n=length(ocEx); 

>> sol=initocmat_AE_AE(m,ocEx{2},1:N+1,iptO,opt); 

>> c=bvpcont(’extremal2ep’,sol,[],opt); 

>> store(m,’extremal2ep’);ocEx=extremalsolution(m);n=length(ocEx); 

The analogous computations are done for the second PCSS 0 (ocEP{4>), yielding a second (heteroge¬ 
neous) indifference threshold point. 
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>> opt=setocoptions(opt,’OCCONTARG’,’MaxStepWidth’,lei, 1 InitStepWidth’,5e-l,’ 
MaxContinuationSteps’,30); 

>> eval=real(eig(ocEP{l}));eval(eval>0)=[];T=10/min(abs(eval)); 

>> sol=initocmat_AE_EP(m, ocEP-fl},1:N+1,ocEP{4}.y(1:N+1),opt,’TruncationTime’ ,T) ; 

>> c=bvpcont(’extremal2ep 1 ,sol,[],opt); 

>> store(m,’extremal2ep’);ocEx=extremalsolution(m);n=length(ocEx); 

>> opt=setocoptions(opt,’OCCONTARG’,’MaxContinuationSteps’,inf); 

>> eval=real(eig(ocEP{17}));eval(eval>0) = [] ;T=10/min(abs(eval)); 

>> sol=initocmat_AE_EP(m,ocEP{17},1:N+1,ocEP{4}.y(1:N+l),opt,’TruncationTime’,T); 

>> c=bvpcont(’extremal2ep’,sol,[],opt); 

>> store(m,’extremal2ep’);ocEx=extremalsolution(m);n=length(ocEx); 

>> elf;xcoord=l;ycoord=l;xvar=’spatialnorm’;yvar=’hamiltonian’; 

>> plotcont(m,xvar,xcoord,yvar,ycoord,’contfield’,’SliceManifold’,’Index’,[6 7] );figure(gcf 

) 

>> opt=setocoptions(opt,’0CC0NTARG’,’InitStepWidth.’,2.5e0,’MaxContinuationSteps’,10); 

>> sol=initocmat_AE_AE(m,ocEx{7},1:N+1,ocEP-fl}.y(1:N+1),opt); 

>> c=bvpcont(’extremal2ep’,sol,[] ,opt); 

>> store(m,’extremal2ep’);ocEx=extremalsolution(m);n=length(ocEx); 

>> ipt=findindifferencepoint(m,6,8); 

>> opt=setocoptions(opt,’0CC0NTARG’,’InitStepWidth’,lei,’MaxContinuationSteps’,inf); 

>> eval=real(eig(ocEP{l}));eval(eval>0)=[];T=10/min(abs(eval)); 

>> sol=initocmat_AE_EP(m,ocEP{l},1:N+1,ipt,opt,’TruncationTime’,T); 

>> c=bvpcont(’extremal2ep’,sol,[],opt); 

>> store(m,’extremal2ep’);ocEx=extremalsolution(m);n=length(ocEx); 

>> sol=initocmat_AE_AE(m,ocEx{7},1:N+1,ipt,opt); 

>> c=bvpcont(’extremal2ep’,sol,[],opt); 

>> store(m,’extremal2ep’);ocEx=extremalsolution(m);n=length(ocEx); 

Finally, we use continuation to find the intermediate indifference threshold points from the transfor¬ 
mation of the first to the second indifference threshold point. 

>> opt=setocoptions(opt,’0CC0NTARG’,’MaxStepWidth’,lei,’InitStepWidth’,5e-l,’SBVPOC’,’ 
BCJacobian’,0); 

>> v=ones(N+l,1); 

>> sol=initocmat_AE_IDS(m,ocEx(9:10),v,ipt0,opt); 

>> c=bvpcont(’indifferencedistribution’,sol,[],opt); 

>> store(m,’indifferencedistribution’); 
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